/*--------------实验二----------------*/
# include <stdio.h>
# include <math.h>

# define PI 3.14

void main()
{
	double Pan_d(double b, double a);//函数声明
	int i, j; //循环变量
	double sum1=0, q_sum2=0, p_sum3=0, san_G[101], san_X[i], san_Z[i], x_k[101], x_DA[180], z_DA[101], Q=2.67 ,G=6.67e-5, M_x=200.0, M_y=200.0, X_yuan=1000.0, Y_yuan=1000.0, N_dian=180.0, R=250.0, k, T;//中间变量
	FILE *fp; //文件指针
	fp = fopen("FuckingAwesome.txt", "w+");
	
	/*------度数化弧度------*/
	 k=(PI/180)*(360/180);   //每一份的弧度
	 
	 /*----坐标计算----*/
	 for(i=0; i<180; i++) {
		 x_DA[i] = R*sin(i*k)+X_yuan;
		 z_DA[i] = R*cos(i*k)+Y_yuan; 	 
	 }
	 /*---核心计算部分-----*/
	 for(i=0; i<101; i++ ) {
		 x_k[i] = (double)i*20;
		 /*-------累加部分---------*/
		 for(j=0; j<180-1; j++) {
			 sum1 = ((1/2*((z_DA[j+1]-z_DA[j])*((x_DA[j]-x_k[i])*z_DA[j+1]-((x_DA[j+1]-x_k[i])*z_DA[j])/pow((z_DA[j+1]-z_DA[j]),2)+pow(((x_DA[j+1]-x_k[i])-(x_DA[j]-x_k[i])),2))*log((pow(x_DA[j]-x_k[i],2)+pow(z_DA[j+1],2))/(pow(x_DA[j]-x_k[i],2)+pow(z_DA[j],2)))) + (((x_DA[j+1]-x_DA[j]-2*x_k[i])*((x_DA[j]-x_k[i])*z_DA[j+1]-(x_DA[j+1]-x_k[i])*z_DA[j])) / (pow(z_DA[j+1]-z_DA[j],2)+pow(x_DA[j+1]-x_DA[j]-2*x_k[j],2)) *(Pan_d(z_DA[j],x_DA[j]-x_k[i]) - Pan_d(z_DA[j+1], x_DA[j+1]-x_k[i])))));
			 q_sum2 = ((( ((z_DA[j+1]-z_DA[j])*(x_DA[j]-x_DA[j+1]-2*x_k[i])) / (pow((z_DA[j+1]-z_DA[j]),2)+pow(x_DA[j]-x_DA[j+1]-2*x_k[i], 2)) ) * (Pan_d(z_DA[j],x_DA[j]-x_k[i]) - Pan_d(z_DA[j+1],x_DA[j+1]-x_k[i]))) - (1/2 * (pow((z_DA[j+1]-z_DA[j]),2)/pow((z_DA[j+1]-z_DA[j]),2)+pow(((x_DA[j+1]-x_k[i])-(x_DA[j]-x_k[i])),2))) * (log((pow(x_DA[j]-x_k[i],2)+pow(z_DA[j+1],2))) /(pow(x_DA[j]-x_k[i],2)+pow(z_DA[j],2))));
			 p_sum3 = ( (pow((z_DA[j+1]-z_DA[j]),2)/(pow((z_DA[j+1]-z_DA[j]),2)+pow(((x_DA[j+1]-x_k[i])-(x_DA[j]-x_k[i])),2))) * (Pan_d(z_DA[j],x_DA[j]-x_k[i]) - Pan_d(z_DA[j+1],x_DA[j+1]-x_k[i])));
			 san_G[i] += (2*Q*G*sum1);
	         san_X[i] += 1/(2*PI) * (M_x*q_sum2 + M_y*p_sum3);
			 san_Z[i] += 1/(2*PI) * (M_x*q_sum2 - M_y*p_sum3);
		 }
		 fprintf(fp, "%lf    %lf    %lf\n", san_G[i], san_X[i], san_Z[i]); 
	 }
	 /*----关闭文件----*/
	 fclose(fp);
}
/*----关于atan(b/a)的判断,并返回一个角度------*/
double Pan_d(double b, double a)
{
	if(fabs(a)>1.0e-15) {
		if(b/a > 0.0) {
			if(a>0.0 && b>0.0) {
				return atan2(b,a);
			}
			else {
				return (-PI + atan2(b,a));
			}
		}
		else {
			return (PI + atan2(b,a));
		}
	}
	else {
		if(b<0.0) {
			return -PI/2;
		}
		else {
			return PI/2;
		}
	}
}